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ABSTRACT 

We calculate the linear momentum flux from merging black holes (BHs) with arbitrary masses 
and spin orientations, using the effective-one-body (EOB) model. This model includes an analytic 
description of the inspiral phase, a short merger, and a superposition of exponentially damped quasi- 
normal ringdown modes of a Kerr BH. By varying the matching point between inspiral and ringdown, 
we can estimate the systematic errors generated with this method. Within these confidence limits, 
we find close agreement with previously reported results from numerical relativity. Using a Monte 
Carlo implementation of the EOB model, we are able to sample a large volume of BH parameter 
space and estimate the distribution of recoil velocities. For a range of mass ratios 1 < m\lmi < 10, 
spin magnitudes of — 0.9, and uniform random spin orientations, we find that a fraction /500 = 
0.12^0 05 of binaries have recoil velocities greater than 500 km/s and /1000 — 0.0271°' 014 have kicks 
greater than 1000 km/s. These velocities likely are capable of ejecting the final BH from its host 
galaxy. Limiting the sample to comparable-mass binaries with mx/ms < 4, the typical kicks are even 
larger, with / 500 = 0.3lt°;^ an d / 1000 = 0.079±g;gg. 

Subject headings: black hole physics - relativity - gravitational waves - galaxies: nuclei 



1. INTRODUCTION 

In the past year there has been remarkable 
progress made in the field of numerical relativ- 
ity (NR). One of the most exciting new results 
is the calculation of the linear momentum flux 
generated by the inspiral, merger, and ringdown 
of black hole (BH) binaries ([Herrmann et al.l 
Baker et al.ll2006l IGonzalez et al.ll2006t 
iKoppitz et al.1 



20061: 



2007T 



2007 



Herrmann et ahl 
.Campanelli et al. 120071 : 
Gonzalez et al .1 120071 : iBaker et al.ll2007f L~ Since the ma- 
jority of this momentum flux is emitted during the 
merger and ringdown, it is difficult to make defini- 
tive predictions for the recoil using only analytic meth- 
ods. However, in the non-spinning case, the post- 
Newtonian (PN) model (jBlanchet et al.l 120051) has pro- 
vided results consistent with NR all along the adia- 
batic inspiral; the effective-one-body (EOB) model can 
reproduce the total recoil, i.e., also the contribution 
from the ring-down phase, but w ith large uncertain- 
ties ([Damour& Gopakumarj |2006|); perturbative mod- 
els (jSopuerta et al.ll2"006D have also spanned the NR pre- 
dictions. The recoil remains an ideal problem for which 
one can benefit enormously from accurate numerical sim- 
ulations. 

Recent NR results predicting very large recoil velocities 
(> 500 km/s) have also called attentio n to the potential 
astrophysical importance o f the recoil (|Campanelli et ail 
[20071: IGonzalez et al.l 12007). For many models of dark 
matter halo growth through hierarchical mergers, super- 
massive BHs at the centers of such halos will inevitably 
also merge unless kicked out of the gravitational po tential 
well from a previous merger (|Menou et al.ll2001l) . The 
NR results have shown that this is indeed possible, but 
have not yet shown whether it is probable or not. 

As an initial investigatio n of this question, we have 
applied the EOB approac h (|Buonanno k Damour|[l999l 
120001: iDamour et "atfeoOOf) to model the BH inspiral and 



ringd own phases including spin-orbit and spin-spin ef- 
fects (jBuonanno et al.1 l2006a| ). described below in § 2. 
Unlike previous analytic approaches, we can now bene- 
fit from a rapidly growing collection of numerical data, 
allowing us to calibrate our model and quantify its uncer- 
tainties. In § 3, we compare a large number of these NR 
simulations to our EOB predictions, agreeing within our 
conservative error estimates. Having established a range 
of confidence in the EOB model, we proceed in § 4 to 
show the results of Monte Carlo simulations with a wide 
range of mass ratios, spin magnitudes, and orientations. 
We thus provide the first estimates of the distribution 
of recoil velocities from BH mergers, and in § 5 briefly 
discuss the astrophysical implications of these results. 

2. ANALYTIC MODEL OF BLACK HOLE RECOIL 

For the two-body dynamics, we use the EOB model 
with spin-orbit and spin-spin terms included through 
1.5PN and 2PN order, respectively, as described in 
iBuonanno et all (|2006ah . The non-spmnmg, conser- 

yative dynami cs are computed through 3PN order 

(|Damour et al.l l2000h and the ra diation-reaction effect s 
are included through 3.5PN order (jBlanchet et al.ll2004h . 
The initial conditions are taken to replicate an adia- 
batic, quasi-circular inspir al beginning at r = 10m . 
Building on preyious wo r ks (IBuonanno" fc Damourll2000l : 
" [2006at |D amour fc GopakumaT l2"006) , 
(2006bj) were able to match the EOB 
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inspiral-plunge waveform to a linear combination of three 
Kerr quasi-normal ringdown (RD) modes, obtaining 
qualitative agreement with the full NR inspiral-merger- 
RD wave. In the model used here, we assume the fi- 
nal B H spin is given by the linear scaling with mass ra- 
tio of IGonzalez et ail (|2006h . After calculating the in- 
spiral and RD dynamics, we determine the linear mo- 
mentum fl ux using t he rad iative multipole moments de- 
scribed in iThorne I <|1980f l . including the leading-order 
radial velocity and spin-orbit contributions to the indi- 
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vidual modes (Kidder 1995). These methods will be de- 
scribed in greater detail in a companion paper. 

The instantaneous transition from inspiral to ringdown 
turns out to be rather sensitive to the point of match- 
ing and it is only partially effective: the gravitational 
wave (GW) frequency resulting from the EOB match- 
ing grows too quickly during the merger-RD transition, 
whereas the NR GW frequency increases more gently. 
As a way of estimating the optimal match point as well 
as the errors associated with this "prompt merger" ap- 
proximation, we calculate the GW energy emitted in 
the inspiral relative to that of the ringdown. Since the 
amplitudes of the individual multipolc modes typically 
increase rapidly with frequency during the inspiral, by 
matching to the ringdown at a later time, the ampli- 
tudes of the excited RD modes and thus the RD energy 
increases significantly. Similarly, if we match too early, 
the resulting RD amplitud es are too small. Following 
iFlanagan fc Hughes! ([1998), we scale the RD energy ac- 
cording to -Erd (/■<<) oc /x 2 /m (we define m — mi + ni2, 
fi = mymilrti) and rj = fx/m, with m\ > m-i) and scale 
the inspiral energy to E- ms {\i) oc jie, where 1— e is the spe- 
cific energy of a test particle at the inner-most stable cir- 
cular orbit (determined here via the effective spin of the 
EOB model). Guided by NR results, we set Ej ns = E RD 
in the equal-mass case (Buonanno et al.ll2006bf ). and use 
the scaling relations to determine the "target" inspiral 
and RD fractions for other BH mass ratios and spins. 

We then vary the matching point in time, calculat- 
ing the recoil for each case, and requiring the frac- 
tion Ei ns /(E[ ns + I?rd) to agree with the target fraction 
within ±15%. This typically corresponds to a range of 
r/m w 2.5 — 3.5 and u) or h ~ 0.1 — 0.2 at the match- 
ing point. For our EOB recoil predictions, we take the 
mean kick velocity over the acceptable match range and 
define l-a errors from the variance over this range. By 
calibrating with the NR results and quantifying its un- 
certainties, the current EOB model can already be used 
to obtain interesting predictions of the recoil velocity. 



3. COMPARISON WITH NUMERICAL SIMULATIONS 

Unlike previous analytic calculations of the recoil, we 
are now in the unique position of being able to compare 
our EOB predictions directly with a growing collection 
of high-resolution NR simulations and thereby establish 
a range of confidence in the results of our analytic model. 
Beginning with the simplest case, in Figure[T]we show in 
the top panel the EOB recoil velocities for non-spinning 
unequal-mass binaries (diamond symbol s with l-a erro r 
bars), along with the classic formula of lFitchett I (Il 983): 
v(r)) ex rj 2 y'l — 477, scaled to the results of lGonzalez et all 
(2006) who find v(?7 ma x) = 175 km/s. In the bottom 
panel, we show the results of equal-mass binaries with 
equal-magnitude spins aligned and anti-aligned with the 
orbital angular momentum. The soli d line is the lin- 
ear sc aling v(a) — 475a km/s given bv iHerrmann et al.l 
( 2007), where we define the dimensionless spin parameter 

Ol,2 = \Sl,2\/ m l,2- 

In Figure [2] we show the combination of the unequal- 
mass and anti-aligned spin effects, with both black holes 
having equal di mensionless spin ma gnitudes a. Following 
the approach of lBaker et al.l (|2007[ ), we plot analytic fits 
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Fig. 1. — Upper: EOB predictions of the recoil velocity for non- 
spinning unequal-mas s BH s, along with l-a errors. The solid line 
follows the IIFitchett 119 83) scaling, normalized to the NR results in 
IGonzalez et alj ^20061 ). Lower: Recoil for equal- mass binaries with 
spins aligned and anti-aligned w ith the orbital angular m omentum. 
The solid line is the linear fit of Herrmann et al. ( 200711 . 



to the EOB predictions of the form 
32V o 2 



v(q,a) 



(i + <z) ; 



^{1-qf + 2(1- q)PK p + K2, (1) 



where q = 1712/1121 and K p = k p a(q + 1). We agree 
closely with their best-fit results, matching Vo = 276 
km/s, and finding slightly higher values of (3 — 1.27, 
and k p = 1.07 (compare to their (3 = 0.84 and k p = 
0.85). This deviation is largest for large spins, where 
they do not have much data and our simplified RD 
matching methods might start to b reak down. F ollow- 
ing the spin-orbit recoil formula in Kidder (1995|), we 
can modify equation ([1]) to include non-planar kicks. 
We write the recoil component in the plane as K p = 
kp(aiCOs9i — qa 2 cos 2 ) and the component out of the 
plane as K z = k z (a\ sin 6\ cos 4>\ — qa,2 sin O2 cos 4>2 ) , with 
#1,2 the angle between the spin and angular momentum 
and 01,2 the azimuthal angle of each spin, measured with 
respect to the binary separation vector. From the Monte 
Carlo simulations (described in the next Section), we find 
k z pa 1.7, but with a large scatter between the EOB pre- 
dictions and the simple analytic model. 

As a final check of our model, in Table Q] we compare 
with t h e more general co n figura tions of Koppit z et alJ 
' 2007jl . ICampanelli et all (l2007f). and IGonzalez et all 
20071) . For the configuration of lGonzalez et alj (|2007l) . 



where the spins are anti-aligned with each other and nor- 
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Fig. 2. — EOB predictions of recoil velocity for a range of 
unequal-mass BHs with spins aligned and anti-aligned with the 
orbital angular momentum. The solid lines are fits to equation 
JTJ. To avoid confusion, we have not plotted error bars, which are 
typically ~ 30 — 40% in this case. 
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Fig. 3. — Cumulative distribution function for kick velocities 
"kick < " f° r equal masses mi = m,2 and random spins. The 
dashed lines represent l-cr confidence limits. 



TABLE 1 

Comparison of NR recoil calculations with EOB model 



mi 
nt-2 


ai 


0,2 


01 


02 


"NR 


"EOB 


Ref. 








(deg) 


(deg) 


(km/s) 


(km/s) 




1 


0.58 


0.0 





180 


100 ± 10 


120 ± 70 


[1] 


1 


0.58 


0.15 





180 


135 ± 20 


160 ± 100 


[1] 


1 


0.58 


0.29 





180 


185 ± 5 


260 ± 80 


[1] 


1 


0.58 


0.44 





180 


215 ± 15 


420 ± 210 


[1] 


1 


0.58 


0.58 





180 


260 ± 15 


340 ± 160 


[1] 


2 


0.8 


0.0 


135 





454 ± 25 


360 ± 150 


[2] 


1 


0.73 


0.73 


90 


90 


2450 ± 250 


1700 ± 400 


[3] 


1 


0.8 


0.8 


90 


90 


2650 ± 300 


1850 ± 450 


[3] 



Refer en ces: [1 ] Koppit z et al.l l|2007fl [2] Campancl li et al.l 
(t2007l') [31 IGonzalez et alj H2007T ) 

mal to the angular momentum vector, the value of the 
final recoil is sensitive to the angle the spins make with 
the orbital velocity vector. We maximize over this angle, 
but are still unable to attain their kick magnitudes of 
~ 2500 km/s. However, we see from these different cases 
that there doesn't seem to be a systematic disagreement 
with the NR results: sometimes the EOB overestimates 
the kick, and sometimes underestimates it, so the ve- 
locity distributions integrated over a wide range of BH 
parameter space should still be reasonably reliable. 

4. RESULTS: RECOIL VELOCITY DISTRIBUTIONS 

Having successfully compared the EOB model with 
a range of NR results, we are now in a position to 
carry out a series of Monte Carlo calculations. The first 
model we consider is that of equal-mass BHs with ran- 
dom spins, uniformly distributed in [0 < a < 0.9] and 
[— 1 < cos 8 < 1]. We construct a probability distribu- 
tion function P(v) for the kick velocities by summing over 
a collection of normal distributions with the mean and 
variance for each individual binary system calculated as 
in §f2] The cumulative distribution function f(v) is given 
by f{v) — P(v')dv' , normalized such that /(0) = 1 
and f(v) is the probability that a random binary in the 
sample will have a recoil larger than v. 

Figure shows this distribution for the equal-mass, 
random spin case for 10 4 binaries, with the dashed lines 
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Fig. 4. — Cumulative probability distribution as a function of 
symmetric mass ratio r). The contours of f(v; if) represent the 
probability of having a recoil velocity greater than v. The spins 
have amplitude a\ = a,2 = 0.9 and random orientation. 



corresponding to l-cr confidence limits. The resulting 
uncertainty is on the order of 50%, which, while admit- 
tedly large, can still provide interesting new constraints 
on astrophysical models of black hole mergers. 

Perhaps a more realistic model samples a uniform 
distribution in mass ratios [1 < mi/1112 < 10] and 
spin orientations, and set a = 0.9 for each black hole, 
based on observational arguments that most supermas- 
sive b lack holes are rapidly spinning (|Yu &: Tremaind 
I2002t lElvis et all 120021 : IWang et all I2OO6I ). For these 
parameters, we find a fraction /500 = 0.12^QQg of bi- 
naries have recoil velocities greater than 500 km/s and 
/1000 = 0.027lg have kicks greater than 1000 km/s. 
However, as we saw in § [3l the largest kicks occur with 
nearly equal masses, so these fractions may be biased 
towards smaller kicks due to the wide range of m\jmi 
considered in this model. To factor out the uncertainty 
in mass distributions, we plot in Figure [4] the velocity 
distribution function as a function of 77. As expected, 
the typical velocities increas significantly for r\ > 0.16. 
Limiting the sample to this range of masses, we find 
/500 = 0.31±g;il and /moo = 0.079±g;g||. 

It is quite possible that black holes in the early uni- 
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Fig. 5. — Upper: Contours of i>5o(»7, a )i defined such that 50% 
of BH binaries with a given mass ratio r\ and spin parameter a are 
expected to have recoil velocities less than vso. Lower: Contours 
of J>9o(»?; o,), defined analogously to vsq. 



verse have actually gained most of their mass by mergers 
as oppose d to accretion, and thus m ay not be rapidly 
spinning (Hug hes fc Blandfordl [2003V ' To explore this 
possibility, we performed another calculation with the 
same range of mass ratios and equal spin magnitudes 
ai = a2, uniformly distributed over [0 < ai^ < 0.9], 
again with random orientations. After calculating the 
resulting distribution function f(v;i], a), we define the 
functions «5o(?7, a) as the velocity below which 50% of 
the predicted kicks lie. Similarly, fgo(7y, a) is greater than 
90% of the kick velocities for a given 77 and a. Contours 
of these functions are shown in Figure \5\ As expected, 
we see that systems with large spins and nearly equal 
masses have the highest kicks. 

We found by inspection that the distribution function 
f(v;r],a) can be well approximated for v < 1500 km/s 



f(v; 77, a) w exp 



-In 10 



on 



(V, a) 



(2) 



Again motivated by maximizing the phenomenological 
equation (fT]), we found a relatively good analytic fit 



V9o(q, a) 



32V q 2 



^1 (1 — q) 2 + k 2 a 2 (q + l) 2 



(3) 



with Vq — 560 km/s and k — 1.3, matching the EOB pre- 
diction of vqq within 25% over the entire range of (77, a) 
investigated, and within 10% for the great majority of 
it. As a check of this method, we applied it to the other 
Monte Carlo samples described above and along with 
equation are able to match closely the integrated dis- 
tribution functions f(v) for v < 1500 km/s. 

5. DISCUSSION 

Along with the recent NR calculations of recoil veloc- 
ities, there has been a great deal of discussion of the 
astrophysical consequences of these kicks. Here we only 
give a brief sum mary of these result s . For a more de- 
tailed review, see I Baker et al.1 (|2006| ) ; ICampanelli et al.1 
(|2007| ). and references therein. Of particular interest 
to our calculation is the observational evidence that 
most of the superma ssive BHs in the low-z universe 
are r apidly spinning (lYu fc Tremaind 120021 : I Elvis et all 
2002; I Wang et al.ll2006ri . and ma y have undergone multi- 
ple mergers with w i\lmi < 3 (lHaehnelt fc Kauffmannl 
l2002t iMerrittJ l2006h . and thus have quite possibly re- 
ceived kicks on the order of 500 or even 1000 km/s. 
This should be enough to eject the final BH from spi- 
ral galaxy bulges an d even most giant elliptical galaxies 
([Merritt et al.ll20"04h [globular clusters may favor more 
extreme mass ra tios, but also have signifi cantly smaller 
escape velocities ([Miller fc H amilton 2002)]. At the same 
time, there is some observational evidence against such 
large kicks (iLibeskind et al.ll"2006h . This apparent dis- 
crepancy may be explained with more detailed models 
for the BH mass distributions in hierarchical merger sce- 
narios, as well as the possibility of using certain preferred 
spin orientations th at would reduce the kick magnitude 
(|Schnittmanl l2004). At the same time, to reduce the 
uncertainties in our predictions, we plan to improve the 
EOB model, notably the details of the matching method 
and the determination of the RD frequencies. 

A.B. acknowledges support from NSF grant PHY- 
0603762 and from the Alfred Sloan Foundation. 
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